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for acoustic wave equation with 
discontinuous coefficients 
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Abstract 


This paper concerns the numerical solution of the acoustic wave equation 
that contains interfaces in the solution domain. To solve the interface prob- 
lems with high accuracy, more attention should be paid to the interfaces. In 
fact, any direct application of a high order finite difference method to these 
problems leads to inaccurate approximate solutions with high oscillations at 
the interfaces. There is however, the possibility of deriving some high order 
methods to resolve this phenomenon at the interfaces. In this paper, a sixth 
order immersed interface method for acoustic wave equation is presented. 
The order of accuracy is also maintained at the discontinuity using the jump 
conditions. Some numerical experiments are included which confirm the or- 
der of accuracy and numerical stability of the presented method. 


Keywords: Interface methods; High order methods; Lax-Wendroff method; 
Discontinuous coefficients; Jump conditions. 


1 Introduction 


In this paper we develop a class of high order numerical methods for wave 
equations with discontinuous coefficients. The class of interface problems in- 
volves many problems of real world applications in Science and Engineering, 
such as Seismology, Ocean acoustics, and Electromagnetic. A nave imple- 
mentation of high order methods fails to achieve high order accuracy and 
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produces some spurious oscillations (Gibbs phenomenon) near the disconti- 
nuity. There are several approaches to deal with such a problem of accuracy 
loss. An efficient method for the simulation of these equations should be 
able to reduce dispersion and dissipation errors in the propagation of the 
solution [16]. Many researchers are interested in high order methods for hy- 
perbolic problems. Long time behavior of the solution of these problems is 
an important challenge to the numerical simulation of such problems. 


A good deal of literature exists on the numerical solution of interface 
problems. Two traditional methods are adding viscosity to the problem and 
using flux limiters [5]. A recent approach is using essentially nonoscillatory 
(ENO) or weighted ENO methods [6]. We distinguish the interface methods 
in two sets with or without use of jump conditions. The first set consists 
of those methods, such as the recent works of Gustafsson and his coworkers 
[3, 4] and Leveque [8] that do not impose jump conditions in the formulation 
of the numerical solution of these equations. These methods are based on 
shock capturing methods and Riemman solvers for conservation laws. The 
second set, whose history goes back to the pioneering work of Peskin [11] for 
the simulation of blood flow in the heart, consists of those methods that use 
some sort of jump condition in their formulation. In fact, to achieve high 
order results, it is recommended to use a special high order method in the 
vicinity of the discontinuity. Derivative matching methods is a related issue 
and have been developed by Driscoll and Fornberg [1] for one and two di- 
mensional Maxwell equations and followed by many others researcher such as 
Zhao and Wei [19], which derived a derivative matching approach based on 
the FDTD schemes for Maxwell equations. This scheme is based on fictitious 
point method and in a vicinity of the interface they introduce original points 
in one side and ghost or fictitious points (unknown values) in the other side 
of the interface. Using the derivative matching conditions the values at the 
ghost points are evaluated. The emphasis in this paper is on the second set of 
methods that use physical jump conditions at the interface which are usually 
easily accessible from the physical properties of the problem. Therefore, we 
perform the numerical solutions to satisfy these jump conditions. The im- 
mersed interface method (IIM) considers a standard method for the regular 
points and imposes a new method for the irregular points to update the solu- 
tion at the next time level with the same accuracy as the standard method. 
The implementation of this new method requires the solution of several small 
linear systems to obtain coefficients of the difference method on the irregular 
points, without imposing any significant computational cost on the calcula- 
tions. The explicit relations between two media (in heterogeneous media) 
through jump conditions eliminate the role of the fictitious points and there- 
for, this method attains the high order results without ghost points. There is 
a related class of methods, known as simplified immersed interface methods, 
that modify explicitly the numerical values at the irregular points [12, 17]. 


The immersed interface method has been discussed for various kinds of 
partial differential equations and its implementation to many real world ap- 
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plications has been successful. A full review of this method for interface 
problems of parabolic and elliptic type is available in the recent book of Li 
and Ito [10]. 

In this paper, a sixth order method for the acoustic wave equation with 
discontinuous coefficients is presented. This higher order method in some 
sense is an improvement on the works of Zhang and LeVeque [18] and Farzi 
and Hosseini [2]. There are several methods for the simulation of the time 
evolution for the wave equation. The Lax-Wendroff method is a simple time 
discretization method for implementation. However, to reduce possible dis- 
persion and dissipation errors one can invoke TVD and WENO methods to 
avoid oscillations near the discontinuity [6, 9]. In this paper, we consider 
a coupled application of Lax-Wendroff and the immersed interface method 
on the interface. The contribution of this paper is given in the next sec- 
tions. The extension of this method to any order is direct and is given in 
Section 4. The stability analysis and implementation of the method for two 
dimensional problems is presented in Section 5. Physical jump conditions are 
demonstrated for the one dimensional acoustic wave equation. 

The explicit and closed form discretization formula is obtained in Section 
2. The theory and numerical results are also developed for piecewise smooth 
coefficients (Sections 3, 4 and 5). The numerical results reported in Section 5 
confirm the efficiency of the method to approximate the solution with a well 
presented behavior of wave propagation and also a high order accuracy at the 
interfaces. The long time behavior of the method is illustrated by Test prob- 
lem 3 in numerical results. The order of the new immersed interface method 
is justified in Test problem 4 with numerical order of accuracy. Numerical 
stability of the method is addressed in Test problem 5. 


2 Acoustic Wave Equation 


Let u(x,t) and p(z,t) be the acoustic velocity and acoustic pressure, re- 
spectively. Then the one-dimensional wave equation can be written as the 
following model problem 


uie=("), aw=(22), () 


and A(a) is a function of position consisting of physical quantities such as 
density p(a), sound speed c(x) and « = pc. We first focus mainly on prob- 
lems in which the density and sound speed are piecewise constant functions 
and have a jump discontinuity at the point x = a, which we call the interface, 


where, 
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ae ={ eee — (3) 


(pre ra S a, 


Then application of our formulation for problems with more general piecewise 
smooth coefficients will be addressed. Typical applications in which this 
assumption is appropriate include long range underwater acoustics, various 
seismological problems as well as electromagnetic problems. Throughout the 
paper we use the symbol At for A(x) with « > a and A~ for A(x) with 
x <a. The same meaning will apply to other matrices. 


In this part, we derive a high order method for this equation and we 
postpone the treatment of the nonsmooth solution at the interface to the 
next sections, where we discuss and derive jump conditions and give an ap- 
proximation that extends the high accuracy of the solution, obtained for the 
smooth regions, to the interface. At the left or right of the point of dis- 
continuity we can use any standard method. The Lax-Wendroff method is a 
simple time evolution explicit method for implementation which uses the val- 
ues at the current time level and does not need any knowledge of previously 
calculated values. This method is based on Taylor series expansion in time 
and substitution of the time derivatives with space derivatives using (1). We 
consider 


U(2;,tn41) ~~ U(x;,tn) pS 
k? OU k4 04U b> OU 

+37 pps (tir tn) ety Boge (tiv tn) + Bl ope (tirtn) (4) 

k8 O°U 

tar aes itn) 


where & is the length of the time steps. Replacing the time derivatives by 
the space derivatives and discretizing the space derivatives we get 


1 1 
U (25, tn41) & U(2;, tn) — kAiQ§U (ej, tn) + sh iQ5 Ua), tn) 


1 

— FB GAQ2U (ej, tn) + GR GQL U (ej, tn) (5) 
1 3 i 

— hh AjQs )U(x;,tn) + aes Sir(a;,tn), 


where A; = A(xz;) and c; = c(z;), (a) is central difference formula of order 


p for 2; [5] and we have used the relation A? = cI, in which I is the 2 x 2 
identity matrix. 


Substitution of gy? as explained in [2], which deals only with the advec- 
tion equation, gives a fully discretized representation of the acoustic equation 
by the following matrix-vector equation 
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U;" =U + ee re (6) 
l=1 


where the T';’s are two by two matrices (coefficient matrices) and can be 
expressed as 


Dj = wA(AA + (4-1), jb =1, 2,...,7, (7) 
with 
w= a A(( — N2¢7)? — N27), 
we = —3A((3 — A?c?)? — 4?c?), 
w3= (A(T — A?e?))?, 
wa = seqArA((6 — A2c?)? — A%c?), 


and w7 = wW1, Wg = W2, Ws = wg. h is the spacial grid step length and \ = z. 


In fact, the high order derivatives are not valid at the interface and con- 
sequently the obtained results are of first order or even less. A detailed 
discussion of this phenomenon of losing the accuracy has already been pre- 
sented by Sei and Symes [14]. So, obviously, on each side of the discontinuity 
of A(x) the method is of order six for the acoustic equation, while at the grid 
points near to the discontinuity the method fails to maintain this order of 
accuracy. 


More precisely, suppose that the interface lies between two adjacent grid 
points with indices J and J+1,ie. 27 <a < x74, then this method works 
well at those grid points at which all the required points to update them 
are located completely on the left or right of the interface, but it fails to be 
accurate at the grid points (irregular points) J — 2,J—1,...,J +3. So, we 
need a new scheme to maintain the same order of accuracy at these irregular 
points. 


For a general piecewise smooth coefficient problem we can derive the same 
method, but in this case the derivatives of the coefficient A(x) will come into 
the difference equation. It should be noted, however, that the coefficient ma- 
trices are not as simple as those appearing in (7), but we can follow the same 
procedure to provide a sixth order method for this case as well. Here with 
an efficient derivation of these formulas we can save more in computations. 
In the following, we present the time derivative approximations in (4) and 
approximations of the corresponding space derivatives 
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(Us)3 nel 
(Urs)? & B! eae +BOQ? UP, 

(Uue)e » COQM UP + COQP Ur + COQP uP, 

(Uses)? © DO ela: + DaQaur + DOQy uP 4 DAQOuUP, 
(Outs)? © E meus + BOQ?) U2 + HOU + BOQMUP 
+ EO QL UF, 

Cnue)t  FOQMUS + FOQ@US + FOQMUS + FOQMUS 
FOQP UPS FOO UN 


where, 


) 

) b 

)=-—BO 4’ — BRA", 
c®@) = BO A- 2B@) A’, 

) 

) 

) 


= _c@ 4! — CR) A" — C8 A”, 

D®@ =-CMA-20@ A’ — 308) A", 

D® =-COA-30@) A’, 

D® =—C®)A, 

BY = apn peat ~ pos = De Ar”. 

E@) = —D A —2D@ A’ — 3D@) A" — 4D Am, 

ES) = —-p@A-—3D®) A’ —6D® A", 

E® = -—p@A-—4D@4 A’, 

E®) = —D@A, 

FQ) =—RO 4’ — BQ A” — BOA” — BAA EO) Am 

F®) =—BO 4 ~— 2B) 4’ — 3B 3) 4” —~ ABA A” — 5EO) AM 

F®) = —FO) A —3E® A’ —6EO A" —10E©) A”, 

FY) =-£S8)A-4E Al —10EO A", 
J=-EYA-5E® A’, 

FO) — F&A, 


A similar standard method has been addressed by Qiu and Shu [13], in which 
the finite difference WENO schemes with Lax-Wendroff time discretization 
for solving nonlinear hyperbolic conservation law systems was developed. It 
uses a WENO scheme instead Qv in (8) and proceeds to the final formula- 
tions. 
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3 Jump Conditions 


We study the jump conditions for a piecewise coefficient problem and then 
consider a general variable coefficient wave equation. A test problem for this 
case is also given in the section on numerical results. 


3.1 Piecewise constant coefficient 


In this section we introduce some jump conditions to serve as a tool for 
developing a new method for the irregular points. By irregular points, we 
mean those grid points that in the process of updating to the next time level, 
use grid points on both sides of the interface. If the interface x = a lies in 
the interval (77,2741) then the irregular points for the given method (6) are 
the grid points J—2,J—1,...,J+3. So, there exist six irregular points and 
the method (6) fails to be accurate at these points. 


For the acoustic wave equation we impose the jump conditions [u] = 0 
and [p] = 0 that can be denoted by a single statement 


[U] =0. (9) 


Using these conditions and the wave equation (1), we obtain the following 
relations at the interface [12], 


OU (at, t) a O*U(a- , t) 


BaF aah , k= 0,1, 252.0 (10) 
where, 
C yar | 10 C xog:| ar 
Do = (>) 01]? Dop+i =(>) 0 ot 
ee 


3.2 A general piecewise smooth coefficient 


For a general piecewise smooth coefficient, by using again the condition [U] = 
0 and imposing the relations [U;] = 0,..., [Usseeee] = 0, we obtain 
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UZ =QU, 
Ue = Su. eo Q3U. eu? 
Of = = QaU, a Qs5U; eas Qe. LxLaD? 
Ons = QU, + QsU, fet QoU. LULU + QioU. LLLX (11) 
Ut = = QuU, + Q12U; se Q13U. LLL ae QiaU. LULL + QisU. LLLLX? 
UF see = Qi1eU, eo Qi7U; jot QisU, LLL a QigU, LLLL + Qao0U, LLLLL 
+Qa1U, LLLLLL 
where 
Qi <> -GA, 
One G?( —~ BYGA-+ Bo), 
Q3 aaa G2 B®) 
Qa = G(-CMQ1 +0 — CPQz), 
2 2 
Q5 = emeren ore e@)). 
Q6 gc? 
Q7 = G4(-DQ, + D® — D@Q, — D?Qv), 
Qs = G4(-D® Q3 — DPQ; + D), 
Qo = G4(-D® Q5 + D®), 
Qi = Gd, 
Qu = G(-BQ, — BQ, — BQ, — BQ, + B), (12) 
Qi2 = G(- £2 Q, — E'Qs — EY Qs + EM), 
Qi3 = G?(-E. BE Q6- EY Ose EO 
Ou = G(—B Qi + E), 
Qi5= Cop, 
Qis = G5(— FY Qi — Fa — FP Qa — FP Or — FP Qu + FY) 
_ 76(_ pA. _ pe) ()y. _ pl) (2) 
Qi7 = G*( Bi: Q3— FY Fi’Q3 — Fe’ Qa t+ Ff ) 
Qis = G9(-FY Qe FO Qs ~ POQs + FO) 
Qio = G9(- Fy Qo — FO'Qu t+ FY) 
Qo = G8(-FO Qi, + F) 
Qu = GS FY) 
where G = (—A*)~! and the other matrices have already been introduced in 


previous sections. 


4 Approximation at the interface for acoustic equation 


We first investigat 


e the approximation of the solution at the irregular points 


for the case of piecewise constant coefficients. At these points we impose the 


same method as ( 


6) and let the coefficients to be determined appropriately. 
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Then we obtain seven unknown 2 x 2 matrices to maintain the sixth order 
accuracy of the method. The details are presented at x, and a similar ar- 
gument is also applied to the other irregular points. The symbol J indicates 
a fixed number corresponding to the interval (#7, 41) that contains the 
interface a. 


Theorem 4.1. If the coefficients of (6) satisfy the following linear system of 
equations 


4 7 
Sealant >) eal Di = FP, (i =0,1,...,6) (13) 
i=1 1=5 
where, 
Tly4 : 
a= (5); =O) do, 6,8 hay es, 7) (14) 
FD = (aya —dA7)' — ah, At Re erence (15) 


then, the method (6) is of order 6 at the irregular point xz. 


Proof. To prove this result we consider the local truncation error at x7 up 
to sixth order 


7 
1 1 1 1 
ds — T, U; fas AU, 7 —kA? LL —K7 ASB Ler a Ae LULL 
Ry a AE We ea 


1 1 


Using the relation A? = c?I we get 
7 


1 1 2 1 22 1 34 
L — k dT iin + (AU, = ghe Use + rg Cc AU per = aa* Cc ere 


1 1 
ke 4A LLLLL wD es LLLLLX Ke 1 
+Ta9 c AU, mg" © Y )z + O(K*) (17) 
Now to proceed with the proof we need to expand each term of (17) up to 
sixth order about x = a. To this end, we distinguish two sets of points in 
first summation 


=: 2, opal ee | _ 1 “ 
Uj; 41-4 =U rU, 5 Ue grt Uwe + pa? Uwe 
dt 3s 1 
—_rpU;, ar eU,, 1<1l<4 18 
T7909"! LLULLXL + 790"! LLULLLX? ann Ue Ss ED ( ) 
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1 1 1 
Uj 4-4 = DU" + r,D\U, or 571 DU ee ot gM DU ree a 5g?) DU eee 
Tis see Hi do. 
—r;D —r/D <Il< 1 
T7990"! bp eanee ot 790"! eu ewes) 5<I< q ( 9) 


where U~ = lim,_,,- U(a,t), and 
Pp 2h pap OO, (LS ct): 


Note that we have used the jump conditions (10) in (19). If we substitute 
(18), (19) and similar expansions for other terms into (17) we obtain L as a 
function of U~, UL, Uz; User: Vives: Vereen 20d Usnerer Therefore, to 
achieve sixth order accuracy we have to force the coefficients of these terms 
to be zero. These systems of matrix equations are exactly the same as (13), 


(14) and (15). 


In theorem 4.1, the unknown 2 x 2 matrices I\7;,/ = 1,2,...,7 will be 
obtained by solving the linear system (13). These linear systems can be 
easily converted to some lower order linear systems. As the matrices Dj, 
j = 1,...,7, are diagonal it is possible to decouple these systems to four 
7 x 7 linear systems; e.g., the first 7 x 7 linear system determines the scalar 
unknowns (['y7),,,/ = 1,2,...,7. In fact, because of this property, there 
are only two different coefficient matrices in these four systems of linear 
equations. These properties are valid at all irregular points. It should be 
mentioned that in the tested numerical problems we did not get any ill- 
conditioning warning due to the coefficient matrices. On the other irregular 
points similar relations can be derived. So, at the grid point J—1 one obtains, 


5 7 
s al yi + S- al y1yDi = Fy, (¢=0,1,...,6) 
I=1 i=6 
where, 
Tl] 
a= (5-1, ((=0,1,...,6, 1=1,2, 31) 


At the grid point J — 2 we have 


6 
S 7 ail y-2 + ai7P y-21Di = Fy, (¢=0,1,...,6) 
1=1 


where, 


an = (F — 2), GSO tiec6 Feat 


At the grid point J +1 we have 
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3 7 
Saal y41Dp' +S > al 411 Sy (i =0,1,...,6) 
l=1 l=4 


where, 


At the grid point J + 2 we obtain 


2 7 


S- ol 7421D;) + S- oul soo = Fy, (i =0,1,...,6) 
= i=3 


where, 


on = (5 +2), Gi 476. GST ta 


At the grid point J +3 we have 


7 


oT 743,1.D;7* + So ail 74341 Se (i =0,1,...,6) 
1=2 


where, 


aun = (5 +3)", (i=0,1,...,6, 7=1,2,...,7). 

The given formulation of the immersed interface method demonstrates the 
possibility of the direct extension of these relations to higher orders. The 
closed form formulas for right hand side matrices (15) are valid for lower 
and higher order formulations. For higher order methods a4 should only be 
replaced with a new one; for example, this element for fourth order method 
is a13. The proof of the following theorem is similar to Theorem 4.1 and so 
we omit it. 


Theorem 4.2. If the coefficient matrices of (6) satisfy the following system 
of matrix equations 


m M-+1 
SCaalyiDp'+ S> oily = Fi, (i =0,1,...,M—1) 
l=1 l=m+1 


where for j < J, F*¥ = F> D7 and forj > J +1, FX = Fj‘. Then, the 
method (6) gives a Mth order approximation of the solution of (1) at irregular 
grid x;. 


Now we extend Theorem 4.1 to the case where the coefficients are piece- 
wise smooth. The local truncation error for a general piecewise smooth co- 
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efficient at x7 can be represented as follows 


7 
=e STE 510 541-4 
l=1 


(TU, + oo +P TOU pe2 + BPO Upeen 


= 


he OU peaee + eT Users) , + O(k®). 


where 

1 1 1 1 1 
TY) = = a oe ae poh soa 
T?) = a ea ae rot ae 
T°) = gots ykD® poh = Te 
pa) — 5D" ee es ” atk Fa, 
T) = me + atk, 

1 

pe) — Fr 


To obtain a sixth order method it is required that the matrices I; satisfy 
the following linear matrix system 


4 7 
os onl + De ail Qo” =Ri, t= 1,2,...,7. 20) 
I=1 I=5 
where 
Q 1,1) = i 
1 1 1 : 
20) _ | 2 30-4 5 
Q Qi + 5iQ2 + GriQa + sar Qr 4 ra riQu + rG7t is, 
1 1 1 
Q@) = Qa4 371s 8s eas war iQ17, 
1 iL 
QQ?) =Qgt+ q@s T 5971 213 + a" /Qis, 
1 1 
Q 5,1) = Qio + 54 + 3971 19: 
1 
Q 61) — Qis + gr 20, 
Q™ = Qar. 


and 
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Ri =0, 

Rg = yT 

Rg = 2vaz4T™ + 2°T), 

Rg = 3vaz4T™ + 6201 4T + 6°T®), 

Rs = 4va3aT + 12v?a024T? + 243.01 4T + 24474), 

Re = 5vag aT + 20v7 04 3kT 4+ 60v2 a4 oT 4+ 12003 a4 17 + 12047), 

Ry = 6va54T™ + 30V? a4 4T® + 1200303 4T + 3600402, 47 
+7208 a1 4T® + 720°T, 


We note that the matrices Q% are diagonal 2 x 2 matrices and it is possible 
to solve the linear systems (20) in a similar way as discussed before. 
At the irregular point 7741 we have 


3 7 
SS aT Qe) + S- ails, = Ri, 1=1,2,...,7. 
1=1 [=4 


where 
Qo) = 7, 
G2 gry srs e ariQ;) + airi@r ae aati! | aiid: 
68) = Qs + 515" arias" aid + ager 
OO) =Qrh4 arias! | wt Qi | Oe, 


7 1 1 

Oe = On ag ariQia i 59" Gis 
re = ‘|: = 

Qe) = Or. = g1t@20 » 

Qo) = Qn - 


Similar formulae can be deduced for other irregular points. 


5 Numerical Results 


In this section, some test problems are given to show the efficiency of the 
derived high order method. The simulation results are given for different 
cases. The test problems Test 1 and Test 2 illustrate general behavior of the 
solution of the acoustic wave equation at the interface. Also, the numerical 
results of the test problem Test 4 ( in two cases 4-1, 4-2) are given for the 
acoustic wave equation with piecewise smooth coefficients. The parameter 
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values, the function f(x), and the CFL number are clearly specified in each 
test problem. The numerical order of accuracy and L; and L.,, errors are 
reported for the test problem Test 3, see Table 1, which verifies numerically 
the long time behavior of this approximation. The test problem Test 5 shows 
the numerical stability of the method. The computational cost for the calcu- 
lation of the coefficient matrices at irregular points is independent of N, the 
number of spatial grid points in the discretization. The coefficient matrices 
can be computed in a couple of milliseconds on a desktop computer. 

The numerical results are given for the following acoustic wave equation 


problem [3]: 
Ifa<a: 
r-a PICL — PrCr r-a 
u(x,t) = t t4 ; 
(x,t) ae (f( 5 ) oa A )) 
r-a PICL — PrCr r—-a 
p(x,t) = f(t ) f(t ) (21) 
Cl picy + PrCr Cl 
and ifa<a: 
9) = 
u(e,t) = ((-=—*), 
PICL T PrCr Cr 
2PrCr rLr-a 
p(z,t) = f(t ), (22) 
PICL T PrCr Cr 


where, f(a) is a smooth function. 


Test 1: In this test, we consider f(x) = e~200(#*-1/2)" and the parameters 
are chosen as a = 0, p; = —1,9, = —1.5, q¢ = 1, c, = 0.5. Figures 1 and 
2 illustrate numerical and exact solutions for u and p with N = 320 and 
CmaxA = 0.8. There are two gaussian pulses going to the right and a 
pulse hits the interface and then transmits with a generated reflecting 
pulse. The CFL number in this case is about one and there is no 
spurious oscillation. 


Test 2: In this test, problem we consider a rather high frequency function 
f(x) = sin(30z) with parameters a = 0, p, = 0.5,p, = 1.0,q = 
0.8,c, = 1.0. Figures 3 and 4 illustrate numerical and exact solutions 
for u and p with N = 320 and cmaxA = 0.8. After hitting the interface 
the magnitudes of velocity and pressure are changed. 


Test 3: To illustrate the long time behavior of this method, we consider 

the jump f(x) = sin(7a), as initial data for acoustic equations with 
a = 5, and the same parameters as in Test 2. The numerical results 
are reported in Table 1 for several values of N at t = 507. The nu- 
merical order of accuracy and L, and DL, errors are presented in this 
Table in which LaxW-IIM denotes the Lax-Wendroff immersed inter- 


face method. 
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... Equation 


u(x,0.4) 


p(x,0.4) 


Figure 1: Test 1: Numerical(.) and exact(-) solutions for w and p at t = 0.4 
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It is well known that a typical high order method provides first order 
results for interface problems[14, 15]. In Table 2, we show the result 
of eliminating of jump conditions and using the high order numerical 
method (6) without interface treatment. We can clearly see that the 
results are at most first order. 


Table 1: (Test 3) The £1, Loo errors and the numerical order of accuracy for LaxW-IIM 
method over the whole interval and the same quantities at irregular points for f(x) = 
sin(7a) at t = 507 are reported 


N LaxW-IIM Irregular points 

fy, error order lL. error order Ly, error order lL. error order 
15 8.37E-001 1.56E-001 2.44E-001 1.06E-001 
30 4.03E-002 4.17 3.32E-003 4.89 6.92E-003 4.89 1.81E-003 5.58 
60 1.11E-003 5.06 4.14E-005 5.46 1.43E-004 5.46 3.05E-005 5.75 
120 3.05E-005 5.12 5.26E-007 5.85 2.37E-006 5.85 4.42E-007 6.04 
240 8.90E-007 5.07 7.18E-009 5.98 3.67E-008 5.98 6.54E-009 6.04 


Test 4: In this test, problem we consider two variable coefficient problems. 


Let us define a general form of the variable coefficient, 


_ (o> filme giz) as a, 
(p(x), c(x)) = { (pt + fo(a), ct + 92(x)), >a. 


(23) 
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u(x,0.8) 


p(x,0.8) 


0 
x 


0.5 


1 


Figure 2: Test 1: Numerical(.) and exact(-) solutions for wu and p at t = 0.8 


(Test 3) The £1, Loo errors and the numerical order of accuracy for method 


in equation (6) over the whole interval and the same quantities at irregular points for 
f(x) = sin(ra) at t = 507, without interface treatment, are reported 


N LaxW-IIM Irregular points 

Ly error order JL. error order Ly error order LL. error order 
15 7.91E+000 2.35E+000 5.17E+000 2.35E+000 
30 =3.84E+000 0.99 8.60E-001 1.57 1.65E+000 1.57 8.16E-001 1.45 
60 4.19EF+000 0.12 3.56E-001 1.14 7.33E-001 1.14 3.56E-001 LAT 
120 4.29E+000 0.03 1.84E-001 0.92 3.84E-001 0.92 1.84E-001 0.94 
240 4.33E+000 0.01 9.24E-002 0.97 1.95E-001 0.97 9.24E-002 0.99 


where f, fe, gi and ge are arbitrary and smooth functions which vanish 


at x =a and p-,pt,c~ 


and c* are constants. 


To illustrate the behavior of the numerical solution near the interface, 
we consider the following two sets of functions and we show the numer- 
ical quality of solution for the first set through some figures and for the 
second set of functions Table 4 is given in which the order of accuracy 
is reported numerically. 


1. The first coefficient set: 
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Figure 3: Test 2: Numerical(.) and exact(-) solutions for u at t = 2.5 


P(x,2.5) 


4 1 1 1 
-1 -0.5 0 0.5 1 
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Figure 4: Test 2: Numerical(.) and exact(-) solutions for p at t = 2.5 
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Table 3: Jumps in the coefficients of Test 4-1 


k | fF 0-) | Mor) | o>) | oO | P07) | OF) | M0) | 0t) 
0 0 0 0 0 0.5 a 0.8 1 
1 1 mak 1 -1 1 -l 1 eal 
2 2 2 2 2 2 2 2 2 
3 -l 1 0 0 -l 1 0 0 
4 -4 -4 0 0 -4 -4 0 0 
5 1 -1 0 0 1 -] 0 0 


The parameters for the first set are the same as in Test 2 and for the 
second one we consider c~ = 1,ct = 14V5 y= = 1, pt = -—6. The 
details of the jump discontinuities of the first coefficient set (24) are 
reported in Table 3. From Table 3 it is clearly seen that the first order 
derivatives of p and c are discontinuous and since c is a polynomial, its 
higher order derivatives become zero while the higher derivatives of p 
are discontinuous. 


Figures 5 and 6 for (24) illustrate the quality of the numerical solutions 
for u and p with N = 320 and cmaxA = 0.8. This behavior confirms 
that the method has been able to successfully capture the solution near 
the interface without any spurious oscillations. 


The numerical order of accuracy and errors for (25), the second coeffi- 
cient set, has been shown in Table 4 for the final time t = 0.5 with the 
following initial data, 


2e—-* —e”, x <0, _ jf 2e-* + e*, a4 <0, 
{cose z>0’ DEO oe x>0. (26) 


It should be mentioned that the results of Table 4 have been obtained 
by implementing our formulations with exact coefficients, confirming 
that the true order of accuracy of the presented method for this type 
of coefficients is also 6. Since the computation of jump conditions for 
the case of piecewise constant coefficients is simple, in practice one 
might prefer to use some approximation of the exact coefficients in the 
implementation. 


Since the obtained approximation is close to the exact coefficient, the order of 
accuracy of the numerical results obtained by the method of this paper should 
be closer to the true 6” order. We have examined the presented method 
using best uniform approximation of degree zero for the coefficients near the 
interface and we obtained an order of accuracy of at least 2. We expect 
to get a solution of higher order of accuracy if the best uniform polynomial 
approximation of a higher degree is used. 
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Figure 5: Test 4-1 
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Figure 6: Test 4-1 
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Table 4: (Test 4-2) The L1, Loo errors and the numerical order of accuracy for LaxW-IIM 
method 
N u p 


I, error order L,. error order Ly error order L,. error order 
16 = 2.45e-006 - 1.35e-006 - 2.66e-005 - 1.35¢e-005 - 
32  2.52e-007 3.28 3.10e-008 5.45 9.21e-007 4.85 2.76e-007 5.61 
64 9.69e-009 4.70 4.71le-010 6.04 2.80e-008 5.04 5.09e-009 5.76 
128 3.26e-010 4.89 7.10e-012 6.05 8.94e-010 4.97 8.74e-011 5.86 
256 1.15e-011 4.82 1.15e-0138 5.94 3.09e-011 4.85 1.49e-012 5.87 


5.1 Numerical stability 


A general stability analysis has no direct answer for the given method. Von 
Neumann stability analysis does not work in this situation, because the co- 
efficients change in a nonsmooth way. So, we are left as an open problem 
verifying stability using either the energy method or the GKS theory [3, 7]. 
But here we illustrate the numerical stability of the method through some 
numerical experiments. The long time behavior of the method, which is very 
important factor in real problems has been illustrated for Test 5. The results, 
also, justify the order of accuracy of the method. In this experiment, we have 
considered several examples of random initial conditions and reported their 
results in Figure 7. In this figure the norm of the solutions u and p are given 
versus 7 showing that they do not increase proportional with e. So, our 
method does not show instability, because a relative growth of the solution 
with respect to this factor is a sign of instability [7]. 

A Von Neumann stability analysis with frozen coefficients provide a rea- 
sonable results on the choice of the Courant number. The influence matrix 
for stability analysis is the following block toeplitz matrix 


Peles Te wie Se 
i PEPE Pe tie. Re 
TI, T3 I+ T4 Ts T¢ Ty 
Tr, TI, Ts i+ T4 Ts T¢ Ty 


Tr, 


ee ek es Be 

T, fo [3 7+T4 
Note that with frozen coefficients, the coefficient matrix A(x) is independent 
of a and therefor the coefficient matrices Ij; are independent of 7 and we 
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eliminate this index for simplicity. For linear stability the eigenvalues of G 
should lie in the unit circle in the complex plane. We numerically locate the 
eigenvalues of G. This matrix depends on A(x) and X. Therefor we report 
the stability results for several values of these parameters. 


Table 5: The norm of eigenvalues of influence matrix for different values of parameters. 
The letter p denotes the periodic boundary conditions 


r c p PCG)  p(Gp) 
1.000 1.000 —1.000 0.543 1.500 
2.000 0.500 —1.500 0.554 0.658 
1.250 0.800 0.500 0.837 1.187 
1.000 1.000 1.000 0.558 1.500 


We remark that the boundary conditions have important role in the sta- 
bility of the problem. In the case of frozen coefficients, i.e. A(x) = AT or 
A(x) = A7, the results are shown in the Table 5. In this case we can choose 
A large enough for different values of A(z). While, a comparison between 
different rows of this table shows that in general the eigenvalues are nonde- 
creasing with variation of parameters. Therefor, for nonsmooth coefficients 
that the situation is more complicated, the inequality max, {|c(x)|,1}# <1 
is a reasonable criteria and numerical tests confirm that this criteria in our 
test problems. 


x 


Test 5: In this test we consider an initial condition u(x;,0) = Rje~® 5)? 
and v(z,0) = 2u(z,0), where R; are uniformly distributed random 
numbers in the interval [0,1]. This example is a variant of a similar one 
dimensional case in [7]. The parameters are c; = 1.0, c, = 0.5, pi = 
2.5, pr = 10.0, N = 1000 and cmaxA = 0.99. The results are given 
in Figure 7, which is a typical test among many other tests. There 
are no noise generation visible near the interface and the norms of the 
solutions do not grow with e. 


5.2 Two dimensional problems 


Implementation of high order interface method for two dimensional acous- 
tic wave equations requires high order jump conditions on the interface. In 
most applications the standard jump conditions are available in the litera- 
ture. Such jump conditions are usually given in the normal and tangential 
directions to the interface. Therefor, we need to define a local coordinate in 
a typical point on the interface to obtain the required approximations at the 
interface(see Figure 8). This is done after transformation of the equation to 
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Figure 7: Test 5: The norm and the solution after a long time 


the new local coordinate system in € — 7 plane. The same formulation in one 
dimensional case will direct us to the set of equations to be solved for the 2D 
and 3D cases. 


6 Conclusions and discussions 


In this paper, we have presented a sixth order immersed interface method for 
acoustic wave equation with discontinuous coefficient. The effect of piecewise 
constant and a more general piecewise smooth coefficients on the derived for- 
mulations has been investigated. We have also provided different numerical 
tests which confirm the efficiency of the method and justify their order of 
accuracy and numerical stability. It should be mentioned that, using jump 
conditions do not impose a considerable computational cost in the calcula- 
tions and one should only solve some low order linear systems to obtain the 
coefficients. In fact, the special treatment of the interface is a preprocessing 
stage in the implementation of immersed interface method and without loss 
of overall speed of computation it is also applicable in the parallel comput- 
ers. In the numerical results, we applied the Lax-Wendroff method for time 
discretization. However, the weighted essentially nonoscillatory(WENO) and 
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Figure 8: Local coordinates in a two dimensional grid with a curve interface 


the total variation diminishing(TVD) methods|6, 9] reduce the possible os- 
cillations in the solution. These methods recently have been added to the 
CLAWPACK software [9] for numerical solution of conservation laws. 


References 


1. Driscoll, T. A. and Fornberg, B. A block pseudospectral method for 
mazwells equations I. one-dimensional case, Journal of Comput. Phys. 
140 (1998) 47-65. 


2. Farzi, J. and Hosseini, S. M. A high order method for the solution of a 
one-way wave equation in heterogeneous media, Far East J. Appl. Math. 
36 (3) (2009) 317-330. 


3. Gustafsson, B. and Wahlund, P. Time compact high order difference 
methods for wave propagation in discontinuous media, SIAM J. Sci. Com- 
put. 26 (2004) 272-293. 


4. Gustafsson, B. and Mossberg, E. High order difference methods for wave 
propagation, SIAM J. Sci. Comput. 26 (2004) 259-271. 


5. Gustafsson, B., Kreiss, H. O. and Oliger, J. Time dependent problems 
and difference methods, first ed., John Wiley and Sons, (1996). 


6. Jiang, G. S. and Shu, C. W. Efficient implementation of weighted ENO 
schemes, J. Comput. Phys., 126 (1996) 202-228. 


7. Larsson, J. and Gustafsson, B. Stability criteria for hybrid difference 
methods, J. comp. phys. 227 (2008) 2886-2898. 


24 


10. 


11. 


12. 


13. 


14. 


15. 


16. 


17. 


18. 


19. 


J. Farzi and S. M. Hosseini 


. LeVeque, R. J. Wave propagation algorithms for multidimensional hyper- 


bolic systems, J. comp. phys. 131 (1997) 327-353. 


. LeVeque, R. J. Finite Volume Methods for Hyperbolic Problems, Cam- 


bridge University Press, Cambridge, 2004. 


Li, Z. and Ito, K. The immersed interface method: numerical solutions 
of PDEs involving interfaces and irregular domains, SIAM, 2006. 


Peskin, C. S. Numerical analysis of blood flow in the heart, J. Comput. 
Phys., 25 (1977) 220-252. 


Piraux, J. and Lombard, B. A new interface method for hyperbolic prob- 
lems with discontinuous coefficients. one-dimensional acoustic example, 


J. Comput. Phys., 168 (2001) 227-248. 


Qiu, J. and Shu, C. W. Finite difference WENO schemes with Laz- 
Wendroff type time discretizations, SIAM J. Sci. Comput. 24 (2003) 2185- 
2198. 


Sei, A. and Symes, W. W. Error analysis of numerical schemes for the 
wave equation in heterogeneous media, Appl. Numer. Math. 15 (1994) 
465-480. 


Symes, W. W. and Vdovina, T. Interface error analysis for numerical 
wave propagation, Comput. Geosci. 13 (2009) 363-371. 


Trefethen, L. N. Group velocity in finite difference schemes, SIAM. Rev. 
24 (1982) 113-136. 


Wiegmann, A. and Bube, K. P. The ezplicit-jump immersed interface 
method: finite difference methods for PDE with piecewise smooth solu- 
tions, SIAM J. Numer. Anal., 37 (2000) 827-862. 


Zhang, C. and LeVeque, R. J. Immersed interface methods for wave equa- 
tions with discontinuous coefficients, Wave Motion, 25 (1997) 237-263. 


Zhao, S. and Wei, G. W. High-order FDTD methods via derivative match- 
ing for Maxwell’s equations with material interfaces, J. Comput. Phys., 
200 (2004) 60-103. 


